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ABSTRACT 

Mass accretion onto (proto-)stars at high accretion rates M* > 10~^ Mq yr"-'^ is expected in massive 
star formation. We study the evolution of massive protostars at such high rates by numerically solving 
the stellar structure equations. In this paper we examine the evolution via disk accretion. We consider 
a limiting case of "cold" disk accretion, whereby most of the stellar photosphere can radiate freely 
with negligible backwarming from the accretion flow, and the accreting material settles onto the star 
with the same specific entropy as the photosphere. We compare our results to the calculated evolution 
via spherically symmetric accretion, the opposite limit, whereby the material accreting onto the star 
contains the entropy produced in the accretion shock front. We examine how different accretion 
geometries affect the evolution of massive protostars. For cold disk accretion at 10"'^ Mq yr~^ the 
radius of a protostar is initially small, J?* ~ a few Rq. After several solar masses have accreted, 
the protostar begins to bloat up and for M^, ~ 10 Mq the stellar radius attains its maximum of 
30 — 400 Rq. The large radius ^100 Rq is also a feature of spherically symmetric accretion at the 
same accreted mass and accretion rate. Hence, expansion to a large radius is a robust feature of 
accreting massive protostars. At later times the protostar eventually begins to contract and reaches 
the Zero- Age Main-Sequence (ZAMS) for Af, ~ 30 M©, independent of the accretion geometry. For 
accretion rates exceeding several 10^"^ Mq yr^^ the protostar never contracts to the ZAMS. The very 
large radius of several 100s Rq results in a low effective temperature and low UV luminosity of the 
protostar. Such bloated protostars could well explain the existence of bright high-mass protostellar 
objects, which lack detectable H II regions. 

Subject headings: accretion - stars: early-type - stars: evolution - stars: formation - stars: pre-main 
sequence 
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1. INTRODUCTION 

Recent studies have revealed that the formation pro- 
cess of massive (> 8 Mq) stars differ in many respects 
from that of low-mass (~ 1 Mq) stars (e.g., Zinnecker 
& Yorke 2007). Although signatures of mass accretion 
are widely observed toward forming massive stars, the 
estimated accretion rates usually exceed 10~^ Mq yr~^, 
which are much higher than the typical value in low-mass 
star formation, 
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where Cs and T is the sound speed and temperature in 
natal molecular cloud cores, and G is the gravity con- 
stant. In the past it has been argued that such high 
accretion rates are required for accretion flow to over- 
come the strong radiation pressure from forming massive 
stars (e.g., Larson & Starrfield 1971; Wolfire & Cassinelli 
1987). However, these results were based on the assump- 
tion of spherically symmetric infall. Nevertheless, there 
are other good reasons to expect high accretion rates 
during the formation of massive stars: Massive stars have 
short Kelvin-Helmholz time scales, once formed they con- 
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sumc the available nuclear fuel quickly, and they will lose 
a significant part of their mass through strong stellar 
winds. Thus, massive stars evolve quickly, even while 
accreting. We can expect the time spent in the main 
accretion phase to be significantly less than the main 
sequence lifetime, presumably of the order of a few dy- 
namical times of the star-forming molecular core. For a 
30 Mq star this provides a lower limit to the average ac- 
cretion rate ~ 10~^ Mq yr~^. Accretion rates may vary 
strongly, however. Thus, phases of even higher accretion 
rates are likely. 

The high accretion rates suggest that the nature of 
massive star formation differs from that of low-mass star 
formation. So far, various formation scenarios, such as 
monolithic collapse of massive dense cores (e.g., Yorke & 
Sonnhalter 2002; McKee & Tan 2002, 2003; Krumholz 
et al. 2007, 2009), and the competitive accretion with 
global infall of the cluster-forming clumps (e.g., Bonnell 
et al. 1998; Bonnell, Vine & Bate 2004; Wang et al. 
2010), have been proposed. The proposed formation sce- 
narios are being assessed with high-resolution observa- 
tions of the massive cores and clumps (e.g., Motte et al. 
2007; Marseille et al. 2008; Zhang et al. 2009; Bentemps 
et al. 2010). 

The evolution of massive protostars at such high rates 
is a critical feature of massive star formation which de- 
termines the effect these stars have on their environment 
(e.g., Zinnecker & Yorke 2007). In our previous work 
(Hosokawa & Omukai 2008; 2009, hereafter Paper I), we 
examined the evolution by numerically solving the inte- 
rior structure of a protostar with a dusty accretion enve- 
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lope under the assumption of spherical symmetry. Our 
calculations show that accreting massive protostars have 
certain characteristic features. At M* = 10"'^ M© yr~^, 
for example, the protostellar radius becomes very large, 
exceeding 100 Rq at maximum. The evolution before 
the arrival to the zero-age main sequence (ZAMS) lasts 
until Af, ~ 30 Mq. At even higher accretion rates 
M* > 3 X 10~^ Mq yr~^, the protostar very abruptly in- 
flates and overtakes material in the inner part the accret- 
ing envelope, thus invalidating the assumption of steady 
mass accretion before reaching the ZAMS. 

The very large stellar radius leads to a low effective 
temperature of the protostar, which in turn explains 
some observational properties of massive protostars in 
Orion KL nebula (Morino et al. 1998; Furuya & Shin- 
naga 2009). 

Observations suggest that mass accretion onto mas- 
sive stars is not spherically symmetric, but rather pro- 
ceeds through massive circumstellar disks (e.g., Cesaroni 
et al. 2007). Signatures of rotating infall are widely 
found toward forming massive stars, sometimes associ- 
ated with massive outflows along the rotation axis (e.g., 
Patel et al. 2005; Beltran et al. 2006a). Recent nu- 
merical simulations also suggest that the massive disk is 
a natural outcome of collapse of a massive dense core 
(Yorke & Sonnhaher 2002; Krumholz et al. 2007, 2009). 
Protostellar evolution via disk accretion has been stud- 
ied for lower-mass protostars with low accretion rates 
Af, < 10-"* Mq yr-i (e.g., Palla & Stabler 1992; Beech 
& Mitalas 1994; Hartman et al. 1997). These authors 
considered the limiting case of "cold" disk accretion, 
whereby most of the stellar surface is able to radiate into 
free space and the accreting material brings only a small 
amount of entropy into the star. They showed that, with 
the low accretion rates, disk accretion slightly reduces 
the stellar radius compared to the spherically symmetric 
cases. 

Maeder and coworkers studied the evolution of ac- 
creting massive protostars in the cold disk limit with 
mass accretion rates increasing with the stellar mass 
up to lO"'' Mq yr"^ (e.g., Norberg & Maeder 2000; 
Behrend & Maeder 2001). Recently, Yorke & Bodcn- 
heimer (2008, hereafter YB08) calculated the evolution 
in the same limit with much higher accretion rates 10"^ 
and 10~^ Mq yr"-'^. Their calculated evolution differs 
from the spherical accretion case in an early phase of 
Af, < 8 Mq. The stellar radius is several 10s Rq with 
spherical accretion, whereas it is a few Rq with disk 
accretion. However, YB08 find that the stellar radius 
rapidly increases soon afterwards and exceeds 100 Rq at 
A^* ~ 10 Mq. The subsequent evolution at Af* > 10 i?© 
is quite similar for both spherical accretion and disk ac- 
cretion. 

In this paper, we examine how the different accretion 
geometries affect the evolution of massive protostars at 
high accretion rates. We explain why the two extreme 
cases lead to the same characteristic features of massive 
protostars, such as the large radius, in spite of the dif- 
ferent evolution in the early phases. To this end, we 
calculate the protostellar evolution in the cold disk ac- 
cretion limit with the numerical code used in Paper I. 
This allows us to extract how the protostellar evolution 
depends on the different accretion geometry independent 



of the effects of different stellar evolution codes. For com- 
pleteness, we also compare the calculated evolution with 
previous results by YB08. 

Finally, we stress that this study is very relevant to 
star formation in the early universe, where high accretion 
rates are also expected. Protostellar evolution with high 
accretion rates has been studied with zero or very low 
metallicities under the assumption of spherically sym- 
metric accretion (e.g., Stabler, PaUa & Salpeter 1986; 
Omukai & Palla 2001, 2003; Hosokawa & Omukai 2009b). 
Recent numerical simulations show that circumstellar ac- 
cretion disks are also expected during star formation in 
the early universe (Yoshida, Omukai & Hernquist 2008; 
Stacy, Grief & Bromm 2010). Tan & McKee (2004) ex- 
amined the evolution of primordial protostars via disk 
accretion using simple analytic models. Here, we calcu- 
late the detailed evolution of a primordial protostar in 
the cold disk accretion limit. 

The organization of this paper is as follows. The 
numerical method and the cases considered are briefly 
described in Section [2] In Section [3l we focus on 
the protostellar evolution at the accretion rate M^ = 
10~^ A^© yr~^, which allows a direct comparison to the 
results of YB08. First, we briefly review the evolution 
via spherical accretion studied in Paper I (Section 13.11) , 
and then consider the evolution via cold disk accretion 
(Section 13. 2|) . Evolution of a primordial protostar is in- 
vestigated in Section 13.31 and subsequently contrasted to 
the present-day cases. Variation of the protostellar evo- 
lution with different accretion rates is studied in Section 
m Finally, Section [5] is devoted to our summary and 
conclusions. 

2. NUMERICAL MODELING 
2.1. Modeling Mass Accretion via Disks 

We examine protostellar evolution by numerically solv- 
ing the stellar structure equations. Case (a) in Figure [T] 
presents the schematic picture of the spherical accretion 
studied in Paper I, where we assumed spherical symme- 
try for an accretion envelope as well as the protostar. 
The numerical method we adopted is summarized as fol- 
lows (e.g.. Stabler, Shu & Taam 1980; also see Appendix 
A of Paper I). For a protostar, we solve the four stellar 
structure equations taking into account mass accretion. 
Free-fall flow is assumed for the optically thin part of 
the accretion envelope. If the accretion flow becomes 
opaque to the stellar radiation before hitting the stellar 
surface, we solve for the structure of this optically thick 
part of the flow together with the stellar structure equa- 
tions of the protostar and the jump conditions at the 
accretion shock front. The accretion shock front is the 
outer boundary of the modeled protostar (shock bound- 
ary). 

For the case of mass accretion via a circumstellar disk 
we consider the limiting case whereby the disk-star con- 
necting layer is geometrically thin and most of the stellar 
surface is unlinked to the accretion flow. Figure [1] depicts 
this situation. In this case, the usual photospheric outer 
boundary condition can be applied for accreting proto- 
stars (e.g., Palla & Stabler 1992), 



2 1 GM, 

-fisf — 77" 



3 Ksf Rl 



(2) 



Evolution of Massive Protostars via Disk Accretion 
(a) spherical accretion (b) cold disk accretion 
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Fig. 1. — Schematic figures of a protostar with different accretion geometry: (a) spherically symmetric accretion, and (b) cold disk 
accretion. In the spherical case (a), the accretion flow directly hits the stellar surface forming an accretion shock front. If the flow is 
optically thick before reaching the surface, the photosphere locates outside the stellar surface. Dust grains in the accretion envelope 
evaporate in a dust destruction front far outside the photosphere. In the cold disk accretion case (b), gas predominantly accretes onto a 
circumstellar disk rather than the star. Accreting material settles onto the stellar surface through a geometrically thin layer (or possibly 
through geometrically thin accretion columns — not shown). Heat brought into the star in the accretion flow radiates freely into space until 
the temperature attains the photospheric value. Most of the stellar surface is unaffected by the accretion flow. The energy radiated away 
by the disk and/or accretion columns before the material settles onto the star (the so-called "accretion luminosity" ) needs to be accounted 
for separately from the intrinsic stellar luminosity. 



iVJ; 



^7TRl<yT^i, 



(3) 



where the suffix "sf indicates the values at the stellar 
surface, and a is Stefan-Boltzman constant. In this limit- 
ing case of disk accretion, accreting materials add to the 
star with the smallest amount of entropy, which is the 
value in the stellar atmosphere unattached to the accre- 
tion flow. Hence, the treatment with equations ^ and 
^ supposes the limit of "cold" disk accretion. Note that 
this cold disk accretion is just a limiting case. Gas accret- 
ing onto the star through the disk is likely to have some- 
what higher entropy especially at high accretion rates 
(Hartmann et al. 1997). The spherical accretion is the 
opposite limit where accreting materials carry the high- 
est amount entropy into the star. Thus, we consider 
two extreme limiting cases with the spherical accretion 
and cold disk accretion. In order to solve the protostel- 
lar structure with cold disk accretion, we modified the 
numerical code used in Paper I to deal with the photo- 
spheric boundary condition. We do not solve the detailed 
structure of the flow connecting the star and disk. 

Different geometries of the accretion flow will lead to 
different structures of the evolving protostar. For exam- 
ple, the average entropy in the stellar interior reflects the 
history of entropy brought into the star with accreting 
material, which in turn depend on the accretion geom- 
etry. For spherical accretion at high accretion rates the 
accretion shock front is embedded inside the stellar pho- 
tosphere and most of the high entropy created at the ac- 
cretion shock front is carried to the stellar surface. This 
entropy is efficiently taken into the stellar interior with 
high accretion rates (Paper I). As a result the proto- 
star has a higher average entropy (compared to a non- 
accreting protostar), which leads to a very large radius 
- exceeding 100 Rq at maximum. With disk accretion, 



however, entropy is mainly generated within the disk by 
viscous dissipation. Before the accreted material reaches 
the stellar surface, a large fraction of the generated en- 
tropy will be transported away by radiation. As a result, 
entropy brought into the star with the accreting material 
is significantly lower than in the spherical accretion case. 
The average entropy within the star s is related to the 
stellar radius as 



i?* oc Af, exp[const. x s]. 



(4) 



This relation is derived by substituting the typical gas 
density and pressure within a star of mass M* and radius 
i?, (e.g.. Cox & Giuh 1968); 
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to the definition of specific entropy of ideal monotonic 
gas. 
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where TZ is the gas constant and /i is the mean molecular 
weight. Equation ([4]) demonstrates that, for the same 
stellar mass, the stellar radius is larger with the higher 
average entropy in the stellar interior. Thus, we can 
naively expect that the stellar radius is reduced with disk 
accretion. Contrary to this expectation, however, calcu- 
lations by YB08, adopting the cold disk accretion, show 
that the radius of a massive protostar exceeds 100 Rq at 
maximum as with the spherical accretion. The goal of 
this paper is to explain why the outcomes are so similar 
among the extreme cases of the spherical accretion and 
cold disk accretion. 
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Table 1. Input Parameters of the cases considered 



Case 


M, (Mo/yr)" 


geometry 


Z= 


M.,0 {Mef 


R,,o (Rer 


reference.'^ 


MD4x3-SDm0.59 


4 X 10"^ 


S 


— ^ 


D 


0.02 


0.3 (0.5) 


32.7 (36.8) 


Section |4] 


MD4x3-S 


4 X 10-3 




s 




0.02 


0.3 


32.7 


Section |4] Paper I 


MD3-D 


10-3 




D 




0.02 


0.1 


3.7 


Sectionl3.2l 


MD3-D-b0.1 


10-3 




D 




0.02 


0.1 


3.5 


Sectionl3.2.2l 


MD3-D-CV 


10-3 




D 




0.02 


1.3 


6.2 


Sectionl3.2.2l 


MD3-S 


10-3 




S 




0.02 


0.05 


15.5 


Section 13. 11 Paper 1 


MD3-S-Z0 


10-3 




S 




0.0 


0.05 


13.0 


Section |3.3| Paper 1 


MD3-SDml 


10-3 


S 


->> 


D 


0.02 


0.05 (1.0) 


15.5 (24.4) 


Section 13.2.21 


MD3-SDm0.1 


10-3 


S 


— 5> 


D 


0.02 


0.05 (0.1) 


15.5 (14.1) 


Section 13.2.21 


MD3-SDm0.1-zO 


10-3 


S 


-!> 


D 


0.0 


0.05 (0.1) 


13.0 (14.1) 


Section I3T31 


MD4-S 


10-4 




S 




0.02 


0.03 


8.8 


Section |4] Paper I 


MD4-SDm0.1 


10-4 


S 


— s> 


D 


0.02 


0.03 (0.1) 


8.8 (5.1) 


Section m 


MD4-D-CV 


10-4 




D 




0.02 


1.3 


6.2 


ApoendixlB.il 


MD5-D-CV 


10-5 




D 




0.02 


1.0 


4.3 


AppendixlB.il 


MD5-S-CV 


10-5 




S 




0.02 


1.0 


4.2 


AppendixlB.il Paper I 



a : mass accretion rate, b : geometry of the accretion flow (D: disk accretion, S: spherical accretion), c ; metallicity, d : mass of initial core 
model, e : radius of initial core model, / : subsections where numerical results of each case are presented, g : the suffix SDinX means the 
geometry is switched from spherical accretion to disk accretion at M» = X Mq. For these cases, stellar mass and radius at the switching 
point are listed in the brackets in the columns of A/*_o and -R,,o. 



2.2. Cases considered 

Table 1 summarizes the cases considered in this study 
and their input parameters. For simphcity the protostel- 
lar evolution is calculated with a constant given accretion 
rate for each case. The adopted accretion rates range 
from 10~^ Mq yr-^ up to 4 x 10"'' Af© yr^^ Evolution 
for the accretion rate M^, = 10""^ Mq yr"^ is studied 
in detail. In addition to evolution with cold disk ac- 
cretion we present several cases with spherical accretion 
for comparison. Evolution via disk accretion at the rates 
Af* = lO"-'^ Mq yr-i and 10"'' Mq yr"! was also studied 
by Palla & Stahler (1992). We consider protostellar evo- 
lution with these low rates to compare our results with 
theirs (see Appendix lB.il) . 

The initial model in each case is constructed following 
Stahler, Shu & Taam (1980) or Palla & Stahler (1991) 
with the adopted accretion rate and boundary conditions 
(also see Appendix A. 2 and B.2 in Paper I). The initial 
stellar mass A/*,o is taken as an arbitrary small value. 
The initial entropy profile in the stellar interior is as- 
sumed as a function of the mass coordinate Af , 



s{M) = s,,o + 13 
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where Sc,o is specific entropy at the stellar center, fee 
is Boltzmann constant, toh is atomic mass unit, and 
/3 (> 0) is a free parameter. We adopt /3 = 1 as a 
fiducial value. Unlike the cases of spherical accretion, 
however, we will see that the protostellar evolution with 
the photospheric boundary changes with different initial 
models even at the same accretion rate (see Section l3.2.2l 
below). We show this by presenting the evolution with 
different initial models with /3 = 0.1 (shallower entropy 
profile, case MD3-D-b0.1), and (3 = (homogeneous en- 
tropy profile, cases with the suffix "cv"). 

The geometry of the accretion flow will vary as a func- 
tion of the stellar mass and radius. Initially, the collapse 
of a gravitationally unstable molecular core is mostly 
spherically symmetric. A hydrostatic object is formed 
when the central regions of the molecular core become 
optically thick. A circumstellar disk is later formed 



as material with increasingly higher angular momentum 
falls toward the newly formed protostar. For these cases 
we start the protostellar evolution calculation assuming 
spherical accretion and later switch to cold disk accretion 
(cases with the suffix "SD"). We also present the evo- 
lution of primordial protostars via disk accretion (case 
MD3-SDm0.1-zO). 

3. EVOLUTION OF MASSIVE PROTOSTARS WITH 
A/, = 10-3 Mq YR-1 

3.1. Spherically Symmetric Accretion 

Here, we briefly review the evolution with spheri- 
cally symmetric accretion studied in Paper I. The up- 
per panel of Figure [2] shows the evolution of the stel- 
lar radius and interior structure at the accretion rate 
M, = 10-3 Mq yr-i (case MD3-S). We deflne the fol- 
lowing four evolutionary phases from the behavior of the 
stellar radius: (I) gradual expansion (AT* < 6 Mq), 
(II) swelling (6 AT© < AT* < 10 Mq), (III) contrac- 
tion (10 Af < Af ,, < 30 Mq), (IV) gradual expansion 
(Af* ^ 30 Mq). a key quantity to understand the va- 
riety of the evolution is the ratio between the accretion 
timescale, 

iacc ^ ^, (8) 

Af ' ^ ' 



and the Kelvin-Helmholtz (KH) timescale, 

GAf2 
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The former is the evolutionary timescale, and the lat- 
ter is the timescale over which a star loses energy by 
radiation. The lower panel of Figure [5] displays the evo- 
lution of these timescales. We see that ^kh significantly 
decreases during phase (II). As a result, the balance of 
these timescales changes from face ^ ^kh in phase (I) to 
face > fKH in phase (III) and (IV). The decrease of fKH 
is caused by a rapid increase of the stellar luminosity L.^,, 
which in turn is a result of the decrease in opacity. In 
most of the stellar interior the opacity's dependence on 
density and temperature is approximated by Kramers' 
law K oc pT~^-^. Thus, as the stellar mass increases, 
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Fig. 2. — Evolution of a protostar via spherical accretion at a rate M, = 10~^ Mq yr~^ (case MD3-S, taken from Paper I). Upper panel: 
Evolution of the interior structure of the protostar. The gray-shaded areas denote convective layers. The hatched areas indicate layers 
with active nuclear burning, where the energy production rate exceeds 10% of the steady rate 0.1LD,st/^* for deuterium burning and 
0.1-L*/M, for hydrogen burning. The thin dotted curves show the locations of mass coordinates M = 0.1, 0.3, 1, 3, 10, and 30 Mq. Lower 
panel: Evolution of the accretion timescale iacc (dashed line), and Kelvin-Helmholz timescale (kh (solid line). In each panel the shaded 
background denotes the four evolutionary phases; (I) adiabatic accretion, (II) swelling, (III) Kelvin-Helmholz contraction, and (IV) main 
sequence accretion phases. 



the temperature in the interior rises and the opacity de- 
creases. 

Figure [3] shows snapshots of radial profiles of entropy 
and luminosity in the early phases (I) and (II). At the 
assumed high accretion rate the entropy is initially ef- 
ficiently brought into the stellar interior with minimal 
radiative loss because of the high opacity during phase 
(I). The instantaneous entropy profile simply traces the 
post-shock values from earlier times (adiabatic accre- 
tion) . The average entropy within the star increases with 
accreted mass, which according to equation (jjj) means 
that the stellar radius should also grow. This is indeed 
the case in phases (I) and (II). As the average density 
decreases and the temperature increases, the opacity de- 
creases; radiative heat transport becomes more efficient 
with increasing stellar mass. 

For M* > 6 Mq (phase II) radiative heat transport 
is efficient enough to modify the entropy distribution 



within the star. The deep interior (where dL/dM > 0) 
loses entropy, whereas the outer surface regions (where 
dL/dM < 0) gain a significant amount of entropy. This 
gain in entropy results in the "bloating up" of the star 
up to > 100 i?0 at maximum. Figure [3] shows that 
the boundary between the heat-losing interior and heat- 
gaining outer layer, namely where dL/dM = 0, moves 
toward the stellar surface with increasing stellar mass. 
Stabler, Palla & Salpeter (1986) called this characteristic 
behavior of luminosity profiles a "luminosity wave" . The 
stellar surface luminosity L* rapidly increases when the 
luminosity wave approaches the stellar surface. After the 
luminosity wave passes through the surface, the star can 
efficiently lose energy by radiation. The star contracts 
to maintain virial equilibrium (Kelvin-Helmholz or KH 
contraction; phase III). The interior temperature rises 
during the contraction. Active hydrogen burning begins 
when the central temperature exceeds 10^ K. After that, 
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Fig. 3. — Radial profiles of the specific entropy and luminosity assuming spherical accretion at a rate M« = 10~ Mq yr~^ (case MD3-S, 
taken from Paper I). Each profile is shown as a function of the mass coordinate M. The two panels display representative profiles during 
the first two evolutionary phases: (I) adiabatic accretion (left), and (II) swelling phase (right). The filled circles on the entropy profiles in 
the left panel indicate the post-shock values. The open circles in the right panel show the values for the bottom edges of the convective 
layers. 



the stellar radius increases following the mass-radius re- 
lation of main sequence stars (phase IV) . 

3.2. Case with Cold Disk Accretion 
3.2.1. Fiducial Model 

We now consider stellar evolution with cold disk accre- 
tion at the same rate M* = 10"^ M© yr"i (case MD3- 
D). First, we discuss the fiducial case with the 0.1 Mq 
initial model adopting /3 — 1 in equation ([T]). The top 
panel of Figure |4] shows the evolution of the stellar radius 
and interior structure. We see that the evolution differs 
from the spherical accretion case for M„ < 10 Mq, as 
prev iously demonstrated by YB08 (also see Appendix 
[R2l for comparison to YB08). For Af, < 5 Mq the 
stellar radius is much smaller than for spherical accre- 
tion. An outer convective zone appears in this phase. 
The protostar then abruptly inflates during the period 
5 Mq < M, < 9 Mq. The maximum radius is ~ 90 Rq 
for M, ~ 10 Mq, which is comparable to the results of 
the spherical accretion case. The stellar radius decreases 
for M* > 10 Mq, and then finally follows the mass-radius 
relationship for main sequence stars when M* > 30 AIq . 
The evolution for Af, ^10 Mq is quite similar to what 
we discussed above for spherical accretion. 

We define the following four phases based on the evo- 
lutionary features: (I) convection (M* < 5 Mq), (II) 
swehing (5 Mq < M, < 9 Mq), (III) KH contraction 
(9 Mq < M, < 30 Mq), and (IV) main sequence accre- 
tion phase (M* > 30 Mq). The top panel of Figure [S] 
shows the evolution of the accretion timescale and KH 
timescale. We note that ^kh significantly decreases at 
M^, ~ 6 Mq, when the protostar rapidly infiates. The 
timescale balance sharply changes from iacc ^ ^kh to 
iacc > ^KH here. As with spherical accretion, we at- 
tribute the sequence of the evolutionary phases from (I) 



to (III) to the inversion of the timescale balance. The 
lower panel of Figure [5] shows that this change is due to 
the rapid increase of stellar luminosity L^, caused by the 
decrease of opacity in the stellar interior with increasing 
mass. The detailed evolution in each phase is explained 
below. 

Convection Phase — The top panels of Figures [2] and |4] 
show that the evolution at M^, < 5 Mq is quite different 
between the spherical accretion case and disk accretion 
case. In the disk accretion case the protostellar radius is 
about one tenth of the value obtained for spherical accre- 
tion. This is explained by equation ^ and the fact that 
the entropy content within the star is much lower with 
disk accretion, a natural consequence of the different ac- 
cretion geometry. With spherically symmetric accretion 
gas settles onto the star with high entropy generated in 
the accretion shock front. For cold disk accretion, on the 
other hand, the entropy of the accreting gas is reduced 
to the value in the stellar atmosphere by the gas' ability 
to radiate into free space. Mass accretion hardly affects 
the average entropy in the stellar interior for the case of 
cold disk accretion. Equation (j4|) shows that the stellar 

— 1/3 

radius decreases according to i?* oc M, in this case. 
We see this decrease for M* < 0.5 Mq in Figure 2] 

Figure 0] shows that deuterium burning (D-burning) 
begins when M* ~ 0.4 Mq. This is in stark contrast to 
the spherical accretion case, where it begins much later 
(M* ~ 6 Mq), as shown in Figure [71 This difference 
reflects the different evolution of the maximum tempera- 
ture within the protostar T^^x (middle panel of Fig. S]) . 
The early rise of T,„ax with disk accretion is explained 
by the relation, 



Lt P M* 
n p R^ 



(10) 
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Fig. 4. — Evolution of a protostar via disk accretion at the rate M« = 10"'^ Mq yr~^ (case MD3-D). Top panel : Evolution of the 
interior structure of a protostar. Illustration of the figure is the same as in the upper panel of Fig. [2| The thin solid line represents 
the evolution of the radius via spherical accretion at the same rate (case MD3-S, taken from Paper I). Middle panel : Evolution of the 
maximum temperature within the star Tmax- Values for spherical accretion at the same rate are plotted with a thin solid line. In this 
figure the shaded background shows the four evolutionary phases; (I) convection, (II) swelling, (III) Kelvin-Helmholz contraction, and (IV) 
main sequence accretion phases. Bottom panel : Evolution of the deuterium concentrations of convective layers /d,cvl/2 (dotted lines), and 
mass-averaged deuterium concentration, /d.av (solid line). 



where T is the typical temperature within the star TZ is 
the gas constant, and fi is the mean molecular weight. 
Equations ([5]) were used again to derive the dependence 
on M^,/R^,. Equation (|10p means that, at the same stellar 
mass, the interior temperature is higher for the smaller 
stellar radius. Since the stellar radius i?* is smaller for 
the disk accretion, the central temperature is higher for 
the same stellar mass. Thus, D-burning begins much 
earlier in the disk accretion case. 

Figure 0] shows that convective zones emerge soon af- 
ter deuterium burning begins. First, a convective zone 
appears near the center for M* ~ 0.5 Mq. However, 
this zone disappears soon, when available deuterium is 
exhausted. The D-burning layer moves to the outer part 
of the star, where fresh deuterium still remains. Another 
convective zone forms there, connecting the D-burning 
layer to the stellar surface, where newly accreted deu- 



terium can be mixed down to the D-burning layer. When 
M* > 1 Mq the total energy production rate by D- 
burning becomes nearly equal to the steady-state rate 
(bottom panel of Fig. [5]), 



Lu,st = MJu = 1500 Lg 



M. 



10-3 Mq/yt 



' [D/H] 
2.5 X 10- 



where Su is the energy available from the deuterium 
burning per unit gas mass. 

For M■^, > 1 Mq, the radiative core extends outward 
and the outer convective layer with its D-burning shell 
becomes geometrically thinner. This evolution is a result 
of the decrease of opacity in the stellar interior; radia- 
tive heat transport becomes more efficient with the lower 
opacity. At the end of phase (I) convection ceases where 
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Fig. 5. — Upper panel: Evolution of the accretion timescale tacc (dashed hne) and Kelvin-Helmholz timescale tKH (solid line). Lower 
panel: Evolution of various contributions to the luminosity of the protostar. The thick solid line depicts the luminosity at the stellar surface 
Lf The total energy production rate by each type of nuclear reaction is shown for 1) deuterium burning (Lji coarse dotted), 2) pp-chain 
(Lpp, dot-dot-dashed), and 3) CNO-cycle (LcNi dot-dashed line). The thin dashed horizontal line indicates the steady deuterium burning 
rate L^ ^t. In each panel the shaded background denotes the four evolutionary phases as in Fig|4] 



the radiative heat transport is efficient enough to carry 
away the heat generated by deuterium burning. Fresh 
deuterium is no longer supphed to the D-burning layer 
which has now become radiative. 

Figure H] shows that the stellar radius gradually in- 
creases after the ignition of deuterium burning. This is 
because deuterium burning enhances the average entropy 
in the stellar interior. The surface convective layer ab- 
sorbs the entropy generated by deuterium burning. Fig- 
ure [n]-(I) clearly shows that the entropy in the surface 
convective layer rises with increasing stellar mass. The 
effect of deuterium burning is also seen in Figure[71 which 
depicts the evolution of the stellar radius when deuterium 
burning is artificially turned off at some point (dashed 
lines). We see that, if deuterium burning is turned off at 
M^, = 3 Mq, the stellar radius ceases to increase. 

Figure [7] also presents the evolution via spherical ac- 
cretion without deuterium burning (dashed line) . We see 
that, unlike the case of disk accretion explained above, 
deuterium burning hardly affects the evolution. A key 
quantity to interpret this difference is the opacity when 
deuterium burning begins. With spherical accretion, 
deuterium burning begins late, when M^ ~ 6 Mq and 
the opacity within the star has already decreased suffi- 
ciently. The radiative heat transport is so efficient that 
the heat generated by deuterium burning is carried away 



before the convection arises. With disk accretion, on the 
other hand, the opacity is still high at the ignition of deu- 
terium burning when M* ~ 0.3 Mq. Deuterium burning 
initiates the convection due to the inefficient radiative 
heat transport. In this case deuterium burning affects 
the stellar evolution until opacity becomes low enough 
owing to increasing the stellar mass. 

Swelling Phase — The stellar radius rapidly increases in 
this phase. However, this is not a homologous expansion 
of the star. Only a small percent of the mass near the 
surface significantly bloats up. For A/* ~ 10 Mq, for 
example, 10% of the total stellar mass fills 90% of the 
radius and 99.9% of the volume. Figure|4]also shows that 
the stellar interior structure changes during the swelling. 
The surface convective zone leaves the D-burning layer 
and moves outward for M* ~ 6 Mq. The radiative core 
covers most of the stellar interior. 

What is the mechanism of the significant swelling? We 
contend that deuterium burning cannot be the primary 
cause, as demonstrated in the top panel of Figure [71 
We see that, even if we completely (artificially) turn 
off deuterium burning, the protostar subsequently in- 
flates in all examined cases. Deuterium burning enhances 
the swelling, however; the earlier deuterium burning is 
turned off, the later the swelling occurs with a smaller 
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Fig. 6. — Radial profiles of the specific entropy and luminosity for disk accretion at the rate M« = 10~^ Mq yr~^ (case MD3-D). 
Illustration of the figure is the same as in Fig. |3] The four panels correspond to the four evolutionary phases, (I) convection (upper left), 
(II) swelling (upper right), (III) Kelvin-Helmholtz contraction (lower left), and (IV) main sequence accretion (lower right) phases. In the 
upper panels the open circles indicate the bottom edges of convective layers. 



maximum radius. 

We examine the mechanism of the swehing by con- 
sidering the evolution after turning off deuterium burn- 
ing. Figure [S] shows snapshots of entropy and luminos- 
ity profiles in each case. For M* < 8 Mq we find the 
characteristic feature of the "luminosity wave" in the lu- 
minosity profiles as in the spherical accretion case (see 
Section \3A\ and Fig. [3]). The heat is transported from 
the interior where dL/dM > to the surface layer where 
dL/dM < 0. The maximum luminosity increases with 
increasing stellar mass. At the same time the heat-losing 
interior extends outward toward the stellar surface. This 
characteristic evolution is caused by the decrease of opac- 



ity within the star. More and more entropy embedded in 
the deep stellar interior is transferred outward by radi- 
ation as the opacity decreases. The surface layer bloats 
up, because it has a high entropy through absorbing part 
of the transported entropy. Figure |8] shows that the lu- 
minosity wave reaches the stellar surface for M* ~ 9 Mq . 
After that, dL/dM > throughout the star and there 
is no entropy-receiving layer. However, the swelling still 
continues with disk accretion. This is because the newly 
accreted gas settles onto the star with the same high en- 
tropy of the stellar atmosphere. As a result, the surface 
layer always has high entropy. On the other hand, the 
stellar interior continues to lose entropy with decreasing 
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Fig. 7. — Effects of deuterium burning on protostellar evolution. The upper and lower panels show the evolution via disk and spherical 
accretion for M, = 10~^ Mq yr~^. In both panels the thick solid line represents the stellar radius. The hatched area denotes a layer with 
active deuterium burning, where the energy production rate is larger than O.ILd st/Af*- In the upper panel the dashed lines show the 
evolution of the stellar radius after turning off deuterium burning at 3, 4, and 6 Mq, respectively, which are indicated by open circles. In 
the lower panel the dashed line shows the evolution with no deuterium burning from the beginning. 



opacity. The significant entropy loss in the stellar inte- 
rior finally overcomes the swelling, and the star begins a 
phase of KH contraction. 

Evolution of luminosity and entropy profiles with deu- 
terium burning is shown in Figure [S] The evolution is 
qualitatively similar to that without deuterium burning 
presented in Figure [5] Deuterium burning increases the 
entropy within the star, which makes the heat transport 
after the opacity decreases more significant. Therefore, 
deuterium burning enhances the swelling. Figure[3]shows 
that with deuterium burning the luminosity wave reaches 
the surface when M* ~ 6 Mq. The stellar surface lumi- 
nosity L* significantly increases at this moment, which 
causes the inversion of the timcscalc balance (Fig. [5]) . 



Later Phases — For M* > 9 Mq the protostar begins to 
contract (tKH < iacc)- The protostar efficiently loses its 
energy by radiation. Figure [S] shows that the entropy in 
the stellar interior decreases during this phase. The max- 
imum temperature within the star soon exceeds 10^ K. 
Hydrogen burning begins first via the pp-chain followed 
by CNO-cycle reactions, which lead to a convective core 
emerges when M^. ~ 20 Mq . The stellar surface luminos- 
ity becomes nearly equal to the total energy production 
rate by CNO-cycle hydrogen burning for Af* c± 30 Mq 
(Fig. [5]). The stellar mass when the protostar reaches 
the ZAMS M^.zAMS is analytically estimated by equat- 
ing the accretion timescale to the KH timescale of ZAMS 
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stars: 



ta 



^KH,ZAMS- 



fl2) 



For 10 Mq < M* < 100 Mq the mass-radius and mass- 
luminosity relationships of ZAMS stars (e.g. Schaller et 
al. 1992) are well fitted by: 



i*,zAMS = 1.4 X 10* Le 



-R*,zAMs = 3.9 Rf. 



M» 



M. 



10 Mq 
0.55 



(13) 



. 10 Mq 

Substituting these relations into equation ([T^ . we get 



(14) 



0.645 



Mh.,ZAMS ^ 30.5 Mq 



10-3 Mq yr- 



(15) 



which agrees well with our calculations. Since the above 
argument is independent of different accretion geome- 
tries, so is M*_ZAMS- 

3.2.2. Dependence on Initial Models 

Next, we examine how the protostellar evolution 
changes with different initial models. To this end we 



consider two initial models that differ from the the fidu- 
cial one. One is a 0.1 Mq star with a shallower initial 
entropy distribution (/3 = 0.1 in equation ([7]); case MD3- 
D-bO.l). The other is a fully convective 1.3 Mq star 
(case MD3-D-cv), which was used by Palla & Stabler 

(1992) to calculate the evolution at M* = lO"'* Mq yr"^ 
via cold disk accretion. Figure [9] shows the evolution of 
stellar radius and maximum temperature in these cases. 
The basic evolution is similar. A protostar undergoes 
swelling and contraction, finally reaching the main se- 
quence accretion phase when M* ~ 30 Mq. However, 
the evolution differs quantitatively for the different ini- 
tial models. For example the maximum stellar radius is 
~ 400 Rq for case MD3-D-b0.1, but ~ 30© for case MD3- 
D-cv (the evolution of the interior structure in these cases 
also differs slightly from that for case MD3-D described 
in Section [3?2l see Appendix 1X1 for details). 

This is in stark contrast to the spherical accretion case, 
where the evolution is almost independent of different ini- 
tial models. For spherically symmetric accretion the evo- 
lution of stellar radius in the adiabatic accretion phase 
is well approximated by the relation (Stabler, Palla & 
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Fig. 9. — Variation of evolution via disk accretion with different initial models with an adopted accretion rate M* = 10~^ Mq yr~^. The 
upper and lower panels display the evolution of the protostellar r adius and maximum temperature within the star. In both panels the solid 
line shows the evolution in the fiducial case explained in Section 13.2.11 (case MD3-D), the dashed line represents the evolution assuming a 
shallower initial entropy distribution (case MD3-D-b0.1), and the dot-dashed line shows the evolution with a fully convective initial model 
(case MD3-cv). For comparison we also display the evolution for spherical accretion at the same accretion rate (case MD3-S) with a thin 
solid line. 



Salpeter 1986): 



i?, ~ 26 Rr. 



Mo 



M* 



0.41 



10-3 M© yr-1 



(16) 



Our calculated evolution with spherical accretion (case 
MD3-S) also obeys this relation before the swelling oc- 
curs. Stabler, Palla & Salpeter (1986) showed that the 
evolution of the radius converges to the above relation 
even beginning with different initial models. With spher- 
ical accretion entropy generated at an accretion shock 
front is taken into the stellar interior along with the ac- 
creted gas. If the initial radius is too large compared 
to equation (|16p . the accretion shock front is weak due 
to the shallow gravitational potential well. The entropy 
produced at the accretion shock front is reduced, which 
decreases average entropy within the star. As a result the 
protostar contracts. If the initial radius is too small, the 
opposite occurs, increasing the radius. With disk accre- 
tion, however, this regulation process does not operate. 
The entropy of the accreted material is set to the value 
in the stellar atmosphere. Since the structure of the stel- 
lar atmosphere differs with different initial models, the 
subsequent evolution also changes (also see Hartmann et 
al. 1997 for low-mass protostars). 
How should we choose the initial model? We infer this 



by considering the evolution of accretion flow toward the 
protostar. Just after a protostar is born, the protostar 
gathers material from its immediate vicinity with low 
angular momentum. The geometry of the accretion flow 
is nearly spherical in this early phase. Afterwards, gas 
originally far from the star falls toward the star with 
increasingly higher angular momentum. A circumstel- 
lar disk forms, and the protostar continues to grow via 
disk accretion. The protostellar evolution in this sce- 
nario is calculated by switching the accretion geometry 
from spherical to disk accretion at some point. 

Figure [10] shows some examples of protostellar evo- 
lution, whereby the accretion geometry is switched for 
Af* = 0.1 Mq and 1 Mq. The protostellar radius in- 
creases in the early spherical accretion phase according 
to equation ([TS]) . The radius subsequently decreases after 
the geometry is switched to cold disk accretion, because 
the accreting gas brings lower entropy into the stellar 
interior. The interior temperature increases as the pro- 
tostar contracts. Active deuterium burning begins when 
the maximum temperature exceeds ~ 10^ K. Figure [TU] 
shows that, the earlier the geometry is switched, the ear- 
lier deuterium burning begins. The subsequent evolution 
is similar to that in our fiducial case MD3-D explained in 
Section r3.2l The maximum radius is ~ 50 Rq in both ex- 
amined cases. The protostar reaches the main sequence 
accretion phase at Af* ~ 30 Mq, which is independent 
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Fig. 10. — Evolution of the stellar radius when the accretion geometry is switched from spherical accretion to disk accretion at a constant 
rate M« = 10~^ Mq yr~^. The thick solid and dashed lines show the evolution when the geometry is changed at M« =0.1 Mq and 1 Mq 
(cases MD3-SDm0.1 and MD3-SDml), respectively. The open circles mark the epoch when the total energy production rate by deuterium 
burning reaches 80% of the steady burning rate Lo.st- The filled circles indicate the epoch when the total energy production rate by 
hydrogen burning exceeds 80% of the luminosity at the stellar surface L« . The thin line show the evolution assuming spherical accretion 
only (case MD3-S). 



of the history of accretion geometry. 

3.3. Comparison with Primordial Protostars 

Rapid mass accretion M* ~ 10"'^ Mq yr"^ is also ex- 
pected during the formation of primordial {Z = 0) stars. 
Here, we compare the evolution of primordial protostars 
to that of present-day {Z — Zq) massive protostars. Fig- 
ure [TT] presents the evolution of the stellar radius with 
metallicitics Z = (case MD3-SDm0.1-zO) and Z = Zq 
(case MD3-SDm0.1). In these cases the accretion geome- 
try is switched from spherical accretion to disk accretion 
at M^, =0.1 Mq. The basic evolution with Z = is the 
same as in the present-day case. After deuterium burn- 
ing begins for M» ~ 0.8 Mq, the protostar undergoes 
the four evolutionary phases of convection, swelling, KH 
contraction, and main sequence accretion. 

However, Figure[TT]also shows some quantitative differ- 
ences between the primordial and present-day cases. At 
Z = 0, for example, the swelling occurs for Af, ~ 6 Mq, 
which is slightly earlier than for solar metalicities. A 
similar difference is also found with spherical accretion, 
as shown in Figure [TTJ These differences are attributed 
to the fact that the opacity is lower at lower metallicity 
in the same thermal state. As explained in Section 13.11 
and 13.21 the swelling occurs as a result of the decrease of 
opacity in the stellar interior with increasing the stellar 
mass. The swelling occurs earlier at Z = 0, because the 
opacity decreases sufficiently at an earlier time. 

Another difference between the primordial and 
present-day case is the stellar radius in the main se- 
quence accretion phase; the ZAMS radius for Z = 
is much smaller than that for Z = Zq, independent of 
the assumed geometry. In the primordial case the initial 
abundance of C, N, and O nuclei, necessary catalysts for 
the CNO-cycle hydrogen burning, is zero. CNO-cycle re- 
actions do not begin until these atoms are provided by 
helium triple-a burning, which occurs at T ~ 10^ K. 



Because pp-chain hydrogen burning does not supply suf- 
ficient energy to halt the star's KH contraction, the pre- 
main contraction phase continues up to the higher central 
densities and temperatures that allow triple-a reactions. 
This ZAMS star thus has a much smaller radius than its 
solar metallicity counterpart. 

4. DEPENDENCE ON MASS ACCRETION RATES 

Finally, we examine how the protostellar evolution 
changes with different accretion rates. As pointed out 
in Section I3.2.2[ the evolution quantitatively varies with 
different initial models even at the same accretion rate. 
In this section we focus on the cases where the accre- 
tion geometry evolves from spherical accretion to disk 
accretion in early phases. Figure [TH displays the evolu- 
tion of the stellar radius at the accretion rates of M* = 
10^-* Mq yr-i (case MD4-SDm0.1), 10"^ Mq yr'^ (case 
MD3-SDm0.1), and 4 x 10"^ ^^ yj.-i (^ase MD4x3- 
SDmO.5). With the rate lO""* Mq jt'^, the protostar 
undergoes the same four evolutionary phases as with the 
rate 10"'^ Mq yr~^. The maximum radius decreases 
with decreasing the accretion rate. The same depen- 
dence is also seen with spherical accretion. Note that, 
however, the maximum radius varies with different ini- 
tial models under the assumption of a freely radiating 
photospheric boundary ("disk accretion") as explained 
in Section [3.2.21 Figure fT2l also shows that, at lower ac- 
cretion rates, the protostar arrives to the main-sequence 
accretion phase at a lower stellar mass. We argue that 
this dependence is robust. The stellar mass at the arrival 
to the ZAMS is just determined by iacc — ^kh.zams as 
shown by equation (|15p , which is independent of different 
accretion geometries and different initial models. 

Figure [T2] shows a qualitatively different evolution at 
the higher rate M* = 4 x 10~^ Mq yr~^ With spher- 
ical accretion a period of rapid expansion occurs when 
Af* ~ 50 Mq. This expansion was originally found by 
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Fig. 11. — Comparison with primordial cases: Evolution of protostellar radius for an accretion rate M = 10~^ Mq yr~^. The thick solid 
and dotted lines show the evolution via disk accretion at Z = 0.02 (case MD3-SDm0.1) and Z = (case MD3-SDm0.1-zO), respectively. 
The accretion geometry was switched from spherical accretion to disk accretion when Af* = 0.1 Mq. The thin solid and dashed lines show 
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Fig. 12. — Evolution of the protostellar radius with disk accretion at different accretion rates Af = 10 ** Mq yr ^ (case MD4-SDm0.1, 
dotted line), 10~^ Mq yr'^ (case MD3-SDm0.1, solid hne), and 4 X 10~^ Mq yr'^ (case MD4x3-SDm0.5, dashed line). The evolution 
with spherical accretion at each rate is also plotted with thin lines for comparison (cases MD4-S, MD3-S, and MD4x3-S). The meanings of 
the open and filled circles on the curves are the same as in Figure [TOl 



Omukai & Palla (2001, 2003) at zero metallicity. The ex- 
pansion occurs when the total luminosity Ltot approaches 
the Eddington luminosity, where Ltot = L* + Lace and 
Lace is an additional component from the accretion shock 
front. Lace = GAf,M»/i?,. In Paper I, we showed that 
the critical Eddington ratio Ltot/^Edd, above which the 
expansion occurs, is about 0.5 at Z = Zq (presented 
again in Figure [T5|) . The critical ratio less than unity is 
because the opacity near the stellar surface is somewhat 
higher than the electron-scattering opacity. 

When Ltot reaches 0.5LEdd in the KH contraction 
phase, rapid expansion occurs and inhibits growth of 



the protostar by steady mass accretion. We can derive 
the critical accretion rate M*_cr, with which the proto- 
star barely reaches the ZAMS without undergoing rapid 
expansion, from the critical Eddington ratio as follows 
(Omukai & Palla 2003; Paper I). In the critical case, 
Ltot reaches O.SLsdd when the protostar arrives to the 
ZAMS. Thus, in the spherical accretion case, the critical 
rate Af* cr satisfies the relation 
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Edd, 
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Fig. 13. — Evolution of Eddington ratio of protostars with disk accretion. The dotted, solid, and dashed lines represent L^/L^^^ for the 
different accretion rates Af* = 10"'* Mq yr'^ (case MD4-SDm0.1), 10~^ M© yr'^ (case MD3-SDm0.1), and 4 X 10~^ Mq yr'^ (case 
MD4x3-SDni0.5), respectively. The thin dashed line shows the evolution of itot/^Edd for spherical accretion with 4 X 10"'^ Mq yr~^ (case 
MD4x3-S) for comparison, whereby Ltot = L* + GMtMt/R*. The thin dotted line indicates the critical ratio L/L^^^ = 0.5. The open 
squares show the luminosity of ZAMS stars Lzams(^*) calculated by Schaller et al. (1992). The filled squares show 3 X LzAMsi^*)- 



where the sufRx "ZAMS" indicates the quantities of 
ZAMS stars. This relation is transformed to 



follows 
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where Kesc is the electron-scattering opacity, a lower limit 
to the actual opacity. The critical rate given by this 
equation is a function of Af*, but its dependence is weak. 
The numerical value is A^*,cr — 3 x 10~^ Mq yr~^ over a 
wide range of Af * , which is consistent with our numerical 
results. 

We see that a similar expansion also occurs for 4 x 
10~^ AfQ yr~^ in the cold disk accretion limit in Figure 
1121 In this case the protostar very abruptly inflates when 
Af* > 20 Mq and never turns to KH contraction. The 
stellar radius finally exceeds 600 Rq for Af* ~ 30 Mq. 
YB08 also reported a similar evolution for the accretion 
rate 10~^ Mq yi~^. Interestingly, it is known that a 
similar expansion also occurs when a main-sequence star 
undergoes rapid mass accretion due to mass exchange in 
close binary systems. Kippenhahn & Mayer-Hofmeister 
(1977) and Neo et al. (1977) considered the accretion 
of deuterium-free gas onto existing main sequence stars. 
The numerical method they adopted is basically the same 
as ours; the protostellar structure was numerically cal- 
culated with a freely radiating photospheric boundary 
condition. 

Figure [13] shows that the critical Eddington ratio for 
the expansion is also i*/fvEdd — 0.5 in the cold disk 
accretion limit. Note that there is no additional lumi- 
nosity component Lace here. The accretion luminosity is 
assumed to have already radiated into free space before 
the accreted material settles onto the star. The criti- 
cal accretion rate in this case is derived as follows. As 
explained in Section [3.21 the stellar luminosity L* signif- 
icantly rises in the swelling phase. The luminosity comes 
from release of gravitational energy after that. Figure [13] 
shows that the stellar luminosity in this phase roughly 



in all cases considered. The stellar luminosity ap- 
proaches izAMs(Af*) just before the arrival to the ZAMS 
at the rate 10"" and 10"^ Mq yr"^ At the rate 
4 X 10~^ Af© yr~^, however, L* approaches 0.5LEdd 
before reaching the ZAMS and the abrupt expansion 
occurs. The stellar mass at which L* ~ 0.5LEdd is 
M* ~ 45 AIq using equation ([T9]) . In order to avoid 
expansion the protostar has to reach the ZAMS with a 
mass Af*_zAMS ^ 45 Mq. This condition is transformed 
to that for accretion rates with equation (1151) . 



M < M* 



2 X 10"^ Mq yr" 



(20) 



The derived critical rate is comparable to that for spheri- 
cally symmetric accretion, as confirmed by our numerical 
calculations. 

Our calculations suggest that a massive protostar sig- 
nificantly bloats up and can not reach the ZAMS by 
steady mass accretion, if the accretion rate is higher than 
a few 10~^ Mq yr^^. This feature is independent of the 
accretion geometry. If mass accretion completely shuts 
off as a result of the fast expansion of the protostar, there 
will be an upper mass limit of pre-main-sequence stars 
around several 10s Mq, as discussed in Paper I. Oth- 
erwise, mass accretion might continue in a non-steady 
fashion. At least, it is certain that the radius of a mas- 
sive protostar reaches several 100s Rq at the high accre- 
tion rate. The very large radius leads to a low effective 
temperature of the protostar. 

For case MD4x3-SDm0.5, for example, the stellar lu- 
minosity exceeds 10^ Lq for Af* > 18 Mq, but the ef- 
fective temperature is only Toff < 10'* K, the value of a 
Af* < 2.5 Mq ZAMS star. The stellar UV luminosity 
is very low for such low Toff- This is important for the 
growth of an H II region around the protostar (e.g., Hoare 
& Franco 2007). Observationally, High-Mass Protostel- 
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lar Objects (HMPOs) are considered to be young forming 
massive stars prior to the formation of an H II region. 
However, a lot of HMPOs have high infrared luminosity 
exceeding 10^ Lq without observable H II regions (e.g., 
Sridharan et al. 2002; Beltran et al. 2006b; Beuther et al. 
2007). There has been some speculation regarding what 
hinders the growth of H II regions around these bright 
sources. ZAMS stars with such high luminosity should 
ionize their surroundings. Some authors have proposed 
that rapid mass accretion quenches the growth of H II 
regions by increasing the gas density and recombination 
rate near the protostar (e.g., Walmsley 1995). With disk 
accretion, however, the H II region would breaks out in 
the polar direction, where the gas density decreases sig- 
nificantly after the material with low angular momentum 
has accreted onto the star or been expelled in an outflow. 
Our calculations well explain the existence of such 
bright HMPOs even with disk accretion. Massive stars 
accreting at high rates bloat up; their UV luminosity is 
too low to form the H II region. 

5. SUMMARY AND CONCLUSIONS 

In this paper we have studied the evolution of mas- 
sive protostars with disk accretion at high rates M, > 
10^^ M© yr~^. We considered the limiting case of "cold" 
disk accretion, whereby most of the stellar surface is not 
affected by the accretion flow, and the accreting mate- 
rial brings a minimum of entropy into the star. We cal- 
culated the evolution of protostars in this limiting case 
by adopting the photospheric (freely radiating) boundary 
condition. The calculated evolution was compared and 
contrasted to the evolution with spherically symmetric 
accretion, corresponding to the opposite limit, whereby 
the accreting gas transports the maximum amount of en- 
tropy into the star. 

First, wc considered in detail the evolution for ikf, — 
10~^ Mq yr~^. The basic evolution via the cold disk 
accretion is summarized as follows: 

The entire evolution is divided into four evolutionary 
phases: (I) convection, (II) swelling, (HI) KH contrac- 
tion, and (IV) the main-sequence accretion phase. The 
evolution in the first three phases varies with different ini- 
tial models. Roughly speaking, the stellar radius reaches 
its maximum 30 — 400 i?© for M» ~ 10 M© at the end of 
the swelling phase (II). On the other hand, the protostar 
begins the final phase (IV) for M, ~ 30 M©, independent 
of different initial models. 

A key physical quantity governing the evolution is the 
opacity in the stellar interior, analogous to the spheri- 
cal accretion case. The opacity decreases with increasing 
stellar mass, which causes the variety of the evolution. 
In phase (I) the opacity is so high that the radiative heat 
transport within the star does not affect the stellar struc- 
ture. In phase (II) the entropy accumulated in the stellar 
interior is gradually transported outward as the opacity 
decreases. A part of the transported entropy is received 
in a thin layer near the stellar surface. Due to the high 



entropy deposited into this layer it significantly bloats 
up, resulting in the swelling seen in phase (II). After- 
wards, the opacity decreases sufficiently, allowing the star 
to lose heat by efficiently radiation. The protostar con- 
tracts and the interior densities and temperatures rise. 
This is the KH contraction phase (HI). The protostar 
finally reaches the main-sequence accretion phase (IV), 
just after the active hydrogen burning begins. 

A major difference between the results for spherical ac- 
cretion and cold disk accretion is the role of deuterium 
burning. In the limit of cold disk accretion the accreting 
gas brings a minimal amount of entropy into the star. 
Since the stellar radius is smaller because of the lower 
average entropy in the stellar interior, the protostellar 
radius is initially small, i?, ~ a few xi?©. However, deu- 
terium burning begins earlier, which raises the entropy 
content within the star. This enhances the swelling. 
Consequently, the maximum radius is as large as that 
with spherical accretion. 

Next, we calculated the protostellar evolution at zero 
mctallicity in the cold disk accretion limit. The pro- 
tostellar evolution varies with metallicities in the same 
manner as with spherical accretion; at Z — 0, the pro- 
tostar swells up slightly earlier, and reaches the ZAMS 
with a smaller stellar radius. 

Finally, we examined the evolution with different ac- 
cretion rates. Our calculations show some dependences 
on accretion rates which were also found with spherical 
accretion. First, the stellar mass at the arrival to the 
ZAMS is higher for higher accretion rates. However, if 
the accretion rate is higher than a few 10~^ M© yr^^, 
the stellar radius continues to increase without return- 
ing to KH contraction. The protostar never reaches the 
ZAMS by steady mass accretion at such high rates. 

The fact that a massive protostar bloats up to a ra- 
dius of several 100s i?© with an accretion rate exceeding 
10""^ M© yr~^ is a very robust result of our studies. The 
large radius leads to a low effective temperature and a 
correspondingly low stellar UV luminosity. Thus, the 
growth of an H II region will be delayed until mass ac- 
cretion onto the star deceases significantly. This explains 
the existence of HMPOs with high luminosities exceeding 

104 Lq. 
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Fig. Al. — Same as Fig|4]but for a different initial model with a shallower entropy distribution (case MD3-D-b0.1). In the upper panel 
the dashed line shows the evolution of stellar radius after deuterium burning is turned off at the point marked by the open circle. 

APPENDIX 
A. EVOLUTION FOR M. = IQ-^ Mq YR-^ WITH SHALLOW INITIAL ENTROPY PROFILE 



In Section [3 .2. 2) we have seen that the protostehar evolution differs with different initial models under the cold disk 
accretion even at the same accretion rate. In this appendix, we add a more detailed explanation of these differences. 
We revisit case MD3-D-b0.1 fSection l3.2.2p . whereby the initial entropy distribution is sha llow er than for the fiducial 
case MD3-D. The principle difference occurs in the first convection phase (I) (see Figure ETj) . A central convective 
core emerges soon after the onset of deuterium burning when M, ~ 0.3 Mq. The convective core grows outward and 
merges with the outer convective layer near the surface: The protostar becomes fully convective for M* ~ 0.4 Mq. 
This is not seen in case MD3-D, where the core remains radiative even after D-burning begins. 

The subsequent evolution also differs. For case MD3-D (see Figure |4|) the radiative core gradually extends outward 
with increasing the stellar mass. The D-burning layer, located at the bottom of the outer convective layer, also moves 
outward, until the protostar swells up (Af* ~ 6 Mq). For case MD3-D-b0.1, on the other hand, deuterium burning 
continues at the stellar center as long as the protostar remains fully convective and freshly accreted gas (containing 
deuteri um) is mixed downward. The deuterium is significantly depleted in this fully convective phase (lower panel 
of Fig. IA1|) . When M* ~ 5 Mq a thin radiative layer (a so-called "radiative barrier") forms within the star, and 
the stellar structure changes significantly. This radiative barrier blocks the convective transport of newly accreted 
deuterium into the stellar interior. Without the refurbishment of deuterium, the central zones quickly consume all 
available deuterium and become radiative. Deuterium burning continues in a thin layer outside the radiative core. 

The evolution described ab ove i s simlar to that of intermediate-mass protostars as calculated by Palla & Stabler 
(1991, 1992; also see Appendix [BT] below). The significant heat input by the deuterium shell-burning triggers the rapid 
swelling of the star. The stellar radius ultimately expands to about 400 Rq. As in case MD3-D, however, deuterium 
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Fig. B1. — Comparison with Palla & Stabler (1992); The evolution of protostcUar radius for our calculations (lines) and Palla & Stahler's 
(symbols). Tbe dashed line and open circles represent the evolution via spherical accretion for 10~^ MQyr~^ (case MD5-S-cv). The 
evolution via disk accretion for 10~^ & 10~* MQyr~^ (cases MD5 & MD4-D-cv) is depicted by the lower & upper solid line and the filled 
circles &; triangles, respectively. The initial models were constructed following Palla & Stabler (1992). 

shell-burning per se is not the main driver of the swelling. Figure ED shows that, even if we turn off D-burning just 
after the swelling begins, the protostar continues to swell up. The protostar expands because a thin surface layer 
bloats up after receiving entropy transported from the inner part of the star. 

B. COMPARISON WITH PREVIOUS WORK 
B.l. Comparison with Palla & Stahler (1992) 

We compare our numerical results with previous calculations by Palla & Stahler (1992), who studied protostellar 
evolution in the limit of cold disk accretion with accretion rates 10~^ Mq yr~^ and 10^'' Mq yr~^. The numerical 
method adopted by Palla & Stahler (1992) is the same as ours (e.g., Stahler, Shu & Taam 1980). The same photospheric 
boundary condition (equations [2] and [3]) was used for the limiting case of cold disk accretion. The initial models they 
adopted are 1.0 M© and 1.3 Mq fully convective stars at accretion rates of 10~^ Mq yr~^ and 10~^ Af© yr~^, 
respectively. We also constructed the same initial models following their method and calculated the subsequent 
evolution at each accretion rate. Figure IbD shows the resulting mass-radius relations. Our calculations agree well with 
the results presented by Palla & Stahler (1992). For these cases the protostar remains fully convective until just before 
the swelling, when a radiative barrier appears as for case MD3-D-b0.1 discussed above. 

In section|4]we discussed another case with 10"'' Mq yr~^ (case MD4-SDm0.1) in addition to that shown in Figure 
IBll Because the initial models differ, the mass-radius relations are slightly different between these cases, although the 
adopted accretion rate and boundary condition are the same. 

B.2. Comparison with Yorke & Bodenheimer (2008) 

Next, we compare the results presented in Section 13.21 to the previous calculations by YB08 (see also discussion 
in Section [T|). YB08 independently calculated the evolution of accreting massive protostars by solving the interior 
structure using a different numerical code than that used by Palla & Stahler (1992). To make the comparison clear 
we recalculate the evolution with their code and briefly analyze the evolution of the stellar interior structure. The 
basic equations are the usual four stellar structure equations, taking into account the effect of mass accretion. The 
adopted boundary condition slightly differs from the one depicted in equations ([2]) and ([3]) , but is essentially the same 
photospheric boundary condition (Bodenheimer et al. 2007; YB08). YB08 solve for a spherical grey stellar atmosphere 
including curvature effects. The stellar interior models are constructed using the Henyey method iterated to match 
the grey atmosphere. The initial model is a fully convective 0.1 Mq star of radius 1 Rq. 

Figure |B2] shows the evolution of the stellar radius and interior structure for M, = 10~^ Mq yr~^. The protostar 
is initially fully convective, until a radiative core appears for Af» ~ 0.8 M©. The subsequent evolution is very similar 
to that for case MD3-D discussed in Section [3.2.11 (see Fig. 2] for comparison). The radiative core gradually grows 
outward as the stellar mass increases. The D-burning layer, located at the bottom of the outer convective layer, also 
moves outward. The stellar radius remains ~ 1 Rq during this phase, which is slightly smaller than for case MD3-D. 
At this smaller radius (and higher central densities and temperatures) central hydrogen burning occurs earlier than 
for case MD3-D. The protostar rapidly swells up when M* ~ 4 Mq. The radiative core simultaneously expands to 
cover most of the star. Central hydrogen burning nearly dies out. 

As noted in Section I3.2.1[ deuterium burning itself is not the cause of the swelling. Figure IB2I shows that even if 
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Fig. B2. — Comparison with Yorke & Bodenheimer (2008): The evolution of protostellar radius and interior structure calculated by 
Yorke & Bodenheimer (2008) for an adopted accretion rate M, = 10"'^ Mq yr~^. The interior structure is illustrated in the same manner 
as in the upper panel of Figure [5] The dashed line represents the evolution of the radius after turning off all nuclear reactions when 
Af* ~ 3.5 Mq, indicated on the plot by the open circle. The dot-dashed line shows the evolution of the radius in case MD3-D (taken from 
Fig.©. 

all nuclear reactions are turned off when A/* ~ 3.5 Mq, the protostar significantly swells up later. After the stellar 
radius reaches its maximum of 120 Rq when M* ~ 13 Mq, the protostar contracts and hydrogen burning increases in 
the center. The convective core also grows due to the increased heat input by hydrogen burning. The star reaches the 
ZAMS for A/* ~ 30 Mq and continues to closely follow the ZAMS as it accretes new material. All these evolutionary 
features are also seen in case MD3-D. We conclude that the evolution presented in Section [3T2l is fully consistent with 
the previous calculation by YB08. 



